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Abnormal diffusion of a single vortex in the two dimensional XY model 
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We study thermal diffusion dynamics of a single vortex in two dimensional XY model. By 
numerical simulations we find an abnormal diffusion such that the mobility decreases with 
time t as f / Int. In addition we construct a one dimensional diffusion-like equation to model 
the dynamics and confirm that it conserves quantitative property of the abnormal diffusion. 
By analyzing the reduced model, we find that the radius of the collectively moving region 
with the vortex core grows as R{t) cx t^^^. This suggests that the mobility of the vortex is 
described by dynamical correlation length as l/\nR{t). 
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1. Introduction 

Vortices, which are topological defects of U(l) symmetry fields, plays an important role in 
low dimensional systems such as thin film super fluids, liquid crystals, layer superconductors 
and Josephson junction arrays(JJAs). A well known example is that two dimensional XY 
(2dXY) model exhibits a vortex driven KT phase transition^ while the elastic theory, which 
ignores the vortices, predicts a unique phase with quasi-long-range order. In this case, i.e., 
a vortices behave as two dimensional Coulomb gas, which make dipole pairs in the ordered 
state. 

Vortex is also an important keyword in dynamical property of the system. In the out-of- 
equilibrium dynamics, relaxation process from certain initial state to the equilibrium at fixed 
temperature environment, it is pointed out that the critical relaxation of the 2dXY model 
seems not to be universal; the dynamical exponent z, which connects dynamical correlation 
length L(t) and time t as L{t) cx depends on the initial state. When the initial state is 

an ordered ground state, z equals 2, which agrees with the result of the elastic description, 
i.e., the Gaussian model. On the other hand z ~ 2.35^ for process quenching from highly 
disordered initial state at high temperature. The latter value of z is, however, considered to 
be a consequence of short time correction caused by topological defects. In the disordered 
initial state the system is filled with vortices and vortex-antivortex pair annihilation is a main 
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process of initial relaxation. This yields logarithmic correction as L{t) (t/lnt)"^/^,^"^ which 
resembles t~i/2-35 -^i^gj^ the observing time is not large enough. The long time asymptotic 
behavior is expressed by z = 2 as well as in the case of the ordered initial state. 

In the phenomena mentioned above, interaction among vortices is important. On the other 
hand, it is also reported that even a single vortex causes abnormal behavior in transport 
property of a JJA system. Under very low magnetic field, which yields very diluted vortices, 
the frequency dependence of the vortex mobility behaves as l/uj.^ This result conflicts with 
the free Coulomb gas picture.®'^ Korshunov explained this experimental observation (deriving 
corresponding facts that the mobility of the vortex decreases as 1/ In t) by analyzing the 2dXY 
model assuming an effective action which involves a memory kernel in the interaction term.® In 
this article, we study the diffusion dynamics of a single vortex based on two models and discuss 
about the origin of the memory effect. At first we show numerical study of the bare 2dXY 
model. Then we derive a one dimensional model as an approximation of the 2dXY model, 
which enables us to understand the phenomena more clearly. Analysis of both models reveals 
that the dynamics of a vortex shows abnormal diffusion, where mean square displacement 
grows as t/lnt. 

2. Numerical analysis on the two dimensional XY model 

We study the XY spin model on a square lattice, whose energy is written as 

S = J^[l-cos(^,-^,)]. (1) 

Here 6i indicates the angle of the XY spin at the i-the site and the summation is taken over 
all nearest neighbor pairs. The dynamics of this model is investigated by the overdamped 
Langevin equation, 

jen.n. 

where <^i{t) is a Gaussian white noise satisfying 

{Ci) = and {Ci{t)Cj{t')) =27]kBT6,,S{t-t'). (3) 

Here (• • • ) means the average over independent noise realizations. In the following, we set the 
coupling constant J, friction coefficient r] and the Boltzmann constant to unity. 

Next let us introduce the quantities to observe. The mean square displacement (MSD) is 
calculated as 

D{t,to)^{m)-Mto)\') (4) 

where X(i) = {X{t),Y(t)) is a position of a vortex at time t and to is a waiting time. The 
velocity auto correlation function is derived from D{t,to) from the relation 

^^L>(t, to) = -2(X(t)X(to)) = -2C{t, to). (5) 
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When the system has time translational symmetry in stationary state, the MSD is a function 
of only t — to and the velocity correlation function is rewritten as 

C{t-to) = l^D{t-to). (6) 

2.1 Simulation settings 

We numerically integrated eq. (2) by the second-order Runge-Kutta method.^ The samples 
are square-shaped including spins and open boundary condition is imposed. In the initial 
state, the phase 9i{t = 0) is given as an angle between the x-axis and the vector — X(0), 
where rj is a position vector of the i-th site and X(0) is the initial position of the vortex core 
set to the center of the sample ((L — l)/2, {L — l)/2). 

Next let us explain how to detect the position of the vortex core. We calculate vorticity 
from the snap shot at time t by summing the phase difference along the closed path of each 
plaquette as 

27rni = [9i+x — 6i] + - ^i+i] (7) 

+ [Gi+y — &i+x+y] + [(^i — &i+y], (8) 

where [x] = x — 27rint(^^) and then [x] returns a value between — tt and tt. Thus can take 
or ±1. A vortex core exists in the plaquette labeled by index i, where nj = 1. The velocity 
of the vortex is calculated as V(t) = (X(f + r) — X(t))/r, where r = mAt is time interval 
between sequential observations and At is an incremental time step of the simulation. We set 
to At = 0.056 and m = 64. When r is set to sufficiently small value, V(t) equals zero in the 
most time steps. This is because the motion of the vortex is discrete and intermittent. In such 
situation we have to note that this velocity averaged for finite time span barely depends on 
the choice of r. 

There are two difficulties in this simulation. One is that pairs of vortex and anti-vortex 
can be created by thermal fluctuation, which makes it very difficult to find the trajectory of 
the vortex initially prepared. Such a thermal excitation, however, is observed with extremely 
small probability at low temperature. We set T = 0.250 ^ 0.28rKT where Tkt is the critical 
temperature of the Kosterlitz-Thouless transition. The other problem is that the vortex feels 
attractive force from the sample edge so that it gets out of the sample in a certain long time. 
In order to observe the long time steady behavior we make an operation to keep the vortex 
around the center of the sample as follows. Let us consider an example case that the vortex 
core moves along the x(right)-direction by a(> 2) steps. At first we delete a columns of spins 
from left edge and then move the all remaining spins to the left by a lattice units. Finally we 
add spins to every empty columns on the right edge by copying the {L — a)-th column in the 
same way. We do similar operation when the vortex core moves to the left, up and down. Such 
operation affects the motion of the vortex core to some extent but this effect will decrease 
with the system size L in a systematic way. 



3/13 



J. Phys. Soc. Jpn. Full Paper 

2.2 Num,erical Result 

If the dynamics of the vortex is the so-called "normal diffusion", the MSD would be 
proportional to Tt and its coefficient means the mobility. In fact D{t) grows slower than 
t- linear behavior as shown in Fig.l. In order to show the deviation from normal diffusion 
apparently we plot Tt/D{t), which can be regarded as a time-dependent effective friction 
coefficient (or inverse of the effective mobility) of the vortex. The results for different waiting 
times to are plotted together. Since the deviation among them is very small as far as t — tQ <^ to, 
the system is considered to be in a stationary state. The coefficient is found to be proportional 
to \nt in long time regime, i.e., 

D{t) (xTt/lnt. (9) 

Finite size effect is observed in the long time limit. 

Figure 1 shows that the coefficients saturate to a finite values which are roughly propor- 
tional to the logarithm of the system size L. 

Next we calculate the velocity auto-correlation function by using the Fourier series of 
the V(t) 

n* — l 

C{nT) = iVfepe-^'^*"'^/^ (10) 
fc=i 

n*-l 

Vfe = ^ V V(to + nr)e2-"'=/^, (11) 

where n* is the total number of observations of the velocity with constant interval r. Although 
the above C{t) barely depends on r, a normalized function C(t)/C(0) does not for t ^ t. 
Note that finite time observation results non- vanishing constant in C{t) for large t (we set 
the observation time equivalent to the waiting time to). Ignoring this finite time effect the 
observation of Fig. 2 supports that 

C(t)oc^-^ = --^(l--^). (12) 
^ ' dt^ Int tlnt2 \ Int J ^ ' 

The correlation function is negative for t > 0. This is natural because the local motion of 

the vortex core driven by random force usually raise the interaction energy with surrounding 

region. Thus restoring force works on the vortex core. (Off course the system has an energy 

invariance against the global translation of spins if boundary effect can be ignored. ) The 

response time should correlate with the range of dragged region. For a single vortex, there is 

no characteristic length scale except the lattice unit and therefore the system has infinitely 

long-time memory. 

3. Reduced model 

In the previous section, we saw that a single vortex does not behave as a Brownian particle 
with normal diffusion but the mobility decreases as l/\nt. For the aim to study the origin 
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Fig. 1. Mean square displacement is plotted with t — t^. The result for three waiting times, five system 
sizes (L=f 6, 32, 64, 128, 256) are shown together. Average is taken over 8192 samples. 



Fig. 2. Vortex velocity auto-correlation function multiplied by t — to- The divergent behavior for large 
t is due to the finite observation time, which is set to be same length with io- 

of the abnormal diffusion the 2dXY model is hard to analyze and numerical simulation is 
rather heavy (note that the observation time of the present simulations is not long enough to 
eliminate the possibility that D{t) oc t^'^ , which is difficult to distinguish from t/\nt for small 



For this reason we propose a reduced model of this system. The fundamental idea is as 
follows. Considering circles with various radii centered on the vortex core, all XY spins in the 
energy minimal state with single vortex are directed the radial direction (see the left figure in 
Fig. 3). We assume that the excited state can be described only the motion of these circles 
and spins on each circle are always along the radial direction (see the right figure in Fig. 3 
and note that spins do not change there position). By taking the center positions as a degrees 




t). 



5/13 



J. Phys. Soc. Jpn. 



Full Paper 




Fig. 3. Schematic diagram of the restricted deformation in the reduced model. We consider virtual 
circles with radius r and center X(r). In fact there are infinite number of circles, whose radius 
varies continuously. The spins direct to the radial direction of the circle on which they are put on. 
The deformation is described only with the positions of the circles. 



of freedom, the resultant equation of motion in continuum limit is 

^X{r, t)=2(-^--^) X{r, t) + Z{r, t). 
dt \ or'' r or J 

{Z{r, t)Z{r', t')) = 2-T5{r - r')5{t - t'). (13) 
vr 

The detail of the derivation is shown in the appendix. Here X{r,t) is the x-component of the 
center of the circle with radius r. The y-component obeys the same equation and decoupled 
with X{r,t). The position of the vortex core is identified with \\m.,f^Q^{r,t). 

4. Numerical integration of reduced model 

To confirm the validity of the above one dimensional model we numerically integrate the 
equation of motion. We write the elastic energy 

2 ^ r,- 2 ^ ^ ' 

j 3 

as a discrete version of eq. (13) where Xi(t) (i = 0, 1, 2, • • • , L — 1) is a degree of freedom 
on the lattice points and AXj = Xj+i — Xi. We use a system with reflective symmetry, i.e., 
rj = j + 1/2 for j < L — 1 and rj = L — j — 1/2 for j > L. Therefore both of Xq and Xl-i 
represents the position of the vortex core. The Langevin equation is written as 



1 



1-1/Aj? 



2i 



+ Zi{t). (15) 



{Zi{t)Zj{t')) = 2Ti5ij5{t - t') (16) 



where 



In this equation temperature can be absorbed by scaling X with yT and therefore we set 
T = 1. We set AX_i = AXj^_i = as an open boundary condition. Actually we approximate 
the denominator in eq. (15), 1 — l/4i^, with unity. 
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Fig. 4. Mean square displacement of the vortex center Xo{t) for the reduced one dimensional model. 




Fig. 5. Vortex velocity auto-correlation function multiplied by i — to for the reduced one dimensional 
model. 



4.1 Abnormal diffusion 

In order to compare the present reduced model with the original two dimensional XY 
model, we calculate the MSD and velocity auto correlation function of Xq^I) and 
Figure 4 and 5 is a result of numerical calculation. The behaviors agree with those of the 
two dimensional XY model; logarithmic correction to the normal diffusion is observed. It can 
be said that the present model holds the essential property of the abnormal diffusion of the 
original model. Furthermore we can perform much longer time simulation on this model than 
on the 2dXY model and observe logarithmic property more clearly. 

4.2 Dynamical correlation length 

On the reduced model we can easily investigate the behavior of off-core region. The MSD 
of i-th site 

Di{t-to) = {\Xi{t)-Xi{to)\^) (17) 
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Fig. 6. Mean square displacement of sites i — 0, 2", 2^, -,2^ are shown together 



is shown in Fig. 6. For small t, Di{t) is proportional to t^^^. This is a standard behavior of an 
stochastic diffusion equation, 



t) = ^X{x, t) + Z{x, t), 



(18) 



which lacks the gradient term in eq. (13). Since local temperature increases with i, the initial 
coefficient of t^/^ term does as well. The core region, however, shows different behavior. Seeing 
in logarithmic scale the growth rate is larger than that for off-core region. This ease to move 
is due to the weak confinement in the vicinity of the free edge on the one side. After Do{t) = 
(Xo(t)^) catches up with Di{t), Di{t) coincides with Do{t). As time goes by, more and more 
regions moves with the core. Therewith the growth Do(^) becomes slower with factor l/\nt 
(seems to become faster in double- logarithmic scale) . This suggests that the dynamics of vortex 
is a collective one and the mobility becomes smaller as its effective radius of the collective 
motion becomes large. 

The range of collective motion can be estimated by the correlation function; 

{Xo{t)X,{t)) 



C{i,t) 



y/Do{t)Di{ty 



A universal scaling function is found so that 

Cii,t) = Fii/Rit))^ 
where 

R{t) oc t^/^. 

In Fig. 7 C{i,t) is plotted as a function of i scaled by t^/"^. 



i/R(t) 



(19) 



(20) 



(21) 



5. Discussions 

We have shown that a single vortex exhibits abnormal diffusion even though there is no 
inter-vortex interaction. What is important is that vortex is not a mere point particle but 
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Fig. 7. The strain correlation function is plotted as a function of site i. The data at different time is 
plotted together, which collapse to a universal curve by scaling horizontal axis by ^/t. 



its motion drags the phase field of the surrounding region. As a natural result velocity auto- 
correlation function has memory with negative correlation. For an isolated vortex, the influence 
of its core motion spreads infinitely large range and then correlation time also diverges. 

Based on the reduced one-dimensional model, clearer picture of the phenomena is obtained. 
In addition, we can perform numerical simulations for large scale both in time and space. We 
found that the correlation length R{t) grows as 

not as (t/lnt)i/2^ Tj^jg 

means that the 

logarithmic correlation does not come from the growth law of the correlation length R{t) itself 
but originates with the size dependence of the mobility as l/lni?(t). This is consistent with 
the finite size effect of the mobility in the 2dXY model. 

The reduced model, eq. (13), has two points that ordinary stochastic diffusion equation 
does not have. Those are the gradient term and the position dependence of temperature. We 
found that the z-linear dependence of the local temperature is not essential for the logarithmic 
correction but the gradient term in eq. (13) is important (not shown here). This term has a 
function to make the deformation amplitude X{r) propagate to the positive direction of r. 
This makes the diffusion of the end of chain slower than that of the ordinary elastic chain. 
This is in contrast with the Bessel equation with n = which has positive sign on the gradient 
term. Since the derivative kernel of the eq. (A-20) is that of the 1st Bessel equation (see the 
appendix), expansion of the solution with the Bessel functions will be useful and the exact 
analysis of the present reduced model may be possible. This is challenging open question. 
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Appendix: Derivation of the reduced model 

A . 1 Elastic energy of single vortex 

At first we consider tfie interaction energy of tlie two dimensional XY model, eq. (1). In 
the elastic continuum approximation, which is justified for the region away from vortex core, 
the energy is written as 

j dr^\\Ve{v)\\ (Al) 

where r = (x^y). A. metastable state having an isolated vortex is written as 

6q{y) = arctan(y/a;). (A-2) 

(Strictly speaking we have to add or subtract tt for x < since the arctangent function returns 
the value between —tt/2 and 7r/2.) The elastic energy of this state is 

Eo = ^J dr2Trri-] =7rln— , (A-3) 

where a is an ultraviolet cut-off length. 

We evaluate the integral in eq. (A-1) supposing 

e{r) = ^o(r') (A.4) 

and r' is related to r as 

r = r' + X(r'), r' = \r'\. (A-5) 

To transform the integral variable from r to r' we first evaluate the Jacobian J to scale 
d\ = JdV: 

d{x',y') \ x'Y' 1 + y'Y'J 

where X' = — -, and r' = (x' = x'/r',y' = y'/r') is the unit vector parallel to r'. To evaluate 
dr 

the gradient 

/dxi dy:_\ 

V%)=^ ^ V'eo(r'), (A-7) 

\ dy dy / 

we need the transform matrix (the Hessian) 

-1 , 



dx' 


dy'\ 


( dx 


dy 


dx 


dx \ 


_ [ dx' 


dx' 


dx' 


dy' 


\ dx 


dy 


dy 


dy/ 


\dy' 


dy' 



J \ -y'X' 1 + x'X' 



(A-8) 



Knowing the gradient with respect to r' 



V'^o(r') = i ( 5 I (A.9) 



x 
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and the integrand of E reads 
dV|V^(r)p -- 



V9{r) 



1 -y'-Y'' 



= dV 
= dV 



Jr' \x' + X' )' 
,l + 2r'-X' + |X'|2 
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(A-IO) 



2 ,1 + r' • X' + |Xf - (r' • X')2 + 0(|X'i3^ 



(A-ll) 



The integration of the first term corresponds to the undistorted vertex energy Eq while that 
of the second term vanishes due to the rotational symmetry with respect to r'. Thus we finally 
find 

E,, = E-E^ = \ Tdr'TrS + odX'p). (A-12) 

Note that two degrees of freedom X{r) and Y[r) are decoupled. Therefore all we have to treat 
is one component field X[r) with one dimensional parameter r. 

By taking variation of the energy function in eq. (A-12) we obtain the energy minimal 
condition 



lE^ 2T,dr-drX{r) = Q. 
oX(r) r 

On the boundary condition; X{G) = Xq and X{R) = 0, the solution is obtained as 



X(r) = Xn 



(A- 13) 



(A.14) 



A. 2 Energy dissipation 

On the next step, we derive effective friction force for X's. By using phase variable the 
energy dissipation rate of the whole system can be written as 

dE^i 



dt 



I dr''\e{r,t)f. 



Remembering tan^(r) = (y — Y) / (x — X), we obtain 

d9{r) sme{r) de{r) cos0(r) 



dx 



dV 



(A-15) 



(A- 16) 



By using this, 



lit 



dX 89 dV 89 
dt dX ^~dtdY 



/ 



dX 



dt 



dr'TT— 





dX 


2 


dY 








+ 




) 




dt 


dt 





sin^ if' + ■ ■ ■ + 0{X^) 
+ 0(X3). 



(A-17) 



Again X and Y are decoupled. Therefore friction force acting on the region (r', r' + dr') is 

dr' (dX dY\ 



r' \dt' dt J 



(A.18) 
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A. 3 Equation of motion 

Wc can construct ovcrdampcd equation of motion for X(r, t) by considering local balance 
between the elastic and friction forces, 

-^X(r,t) = 2Trdr-drX(r,t). 
r at r 

'^Xir,t) = 2(^-l^]xir,t). (A.19) 



dt ' \dr'^ r dr ^ 

By putting u(r, t) = X{r, t)/r, we obtain 

l<'(M)=(^ + ^|-i)»(M). (A.20) 
Assuming separable solution u{r,t) = T{kr) exp{—kH), we obtain the Bessel equation 

with n = 1. It is natural that deformation is not isotropic (n = 0) but anisotropic (ra = 1) 
since the translational motion of vortex breaks the circular symmetry. 
A. 4 Thermal noise 

Finally we consider the property of thermal noise. We introduce Gaussian white noise 

This force Z{r,t) have to satisfy the fluctuation-dissipation relation 

{Z{r, t)Z{r', t')) = 2-TS{r - r')S{t - t') (A-23) 

TT 

to realize the canonical distribution at temperature T in equilibrium. The deviation of this 
random force is not uniform in space but proportional to r. It can be said that effective 
temperature becomes higher with r. 
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